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Abstract 

In  this  report  we  propose  a  state  space  model  of  a  direct  methanol  fuel  cell.  This  dynamical  model  allows  analysis  and  optimisation  of  the  fuel 
cell  in  a  control-theoretical  framework,  for  instance  to  improve  the  overall  dynamics  of  the  fuel  cell  or  to  reject  environmental  disturbances. 
The  model  is  particularly  well  suited  to  synthesise  controllers,  which  are  necessary  to  ensure  stable,  robust  and  efficient  operation  of  a  fuel 
cell.  To  yield  a  state  space  model  of  the  direct  methanol  fuel  cell,  a  system  of  nonlinear  partial  differential  equations  is  transformed  into  a 
system  of  ordinary  initial  value  problems  by  Laplace  transforms,  approximation  in  function  spaces  and  the  method  of  weighted  residuals. 

©  2005  Elsevier  B  .V.  All  rights  reserved. 
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1.  Introduction 

As  a  result  of  the  increasing  usage  of  accumulators,  explor¬ 
ing  alternatives  to  conventional  batteries  is  becoming  more 
and  more  interesting.  A  promising  possibility  is  the  concept 
of  fuel  cells.  Fuel  cells  are  expected  to  reach  higher  degrees 
of  efficiency  because  of  the  direct  conversion  of  chemical 
energy  into  electrical  energy.  Furthermore,  they  produce 
little  noise  and  do  not  need  large  down  times  caused  by 
charging. 

The  direct  methanol  fuel  cell  (DMFC)  is  especially 
promising  for  mobile  applications  with  strongly  limited 
space  and  weight.  Compared  to  other  fuel  cells,  it  has  a 
simple  system  design  which  leads  to  little  required  space. 
Currently,  some  applications  in  laptops,  sailing  boats  etc. 
have  already  been  reported.  Yet  the  existing  DMFCs  are  not 
competitive  with  respect  to  cost/performance  ratio.  In  order 
to  optimise  existing  cells  it  is  important  to  analyse  and,  if 
necessary,  to  control  its  dynamical  behaviour,  e.g.  [1,8]. 

While  general  models  of  a  DMFC  have  already  been 
formulated,  many  of  them  solely  describe  the  steady  state 
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behaviour.  At  the  Institute  of  Chemical  Process  Engineering 
at  Clausthal  Technical  University,  Germany,  a  dynamical 
model  for  a  DMFC  has  been  set  up.  This  model  describes 
methanol,  oxygen,  carbon  dioxide,  temperature,  ionic  and 
electrical  potential  in  the  fuel  cell  and  allows  dynamical 
simulations  with  variable  temperature  for  one-  and  two- 
dimensional  problems  [7]  by  means  of  partial  differential 
equations  (PDEs).  Based  on  these  PDEs,  a  finite  element 
model  (FEM)  has  been  developed  to  analyse  and  optimise 
a  particular  DMFC  operated  at  the  same  institute.  As  a 
result,  it  is  possible  to  simulate  the  reaction  of  the  fuel  cell 
to  certain  changing  conditions,  particularly  to  the  anodic 
concentration  of  methanol  and  to  the  overall  voltage. 

In  order  to  optimise  the  DMFC  it  is  necessary  to  improve 
its  dynamics.  Control  theory  is  the  approved  framework  to 
deal  with  the  dynamics  of  self-acting,  aimed  manipulation 
of  the  inputs  to  achieve  a  desired  output.  Control  theory  also 
provides  heaps  of  standard  procedures  to  design  a  controller 
for  a  certain  input/output  behaviour.  Nevertheless,  it  is  nec¬ 
essary  to  transform  the  partial  differential  equations  of  the 
mathematical  model  into  an  equivalent  standard  form,  the 
state  space  model,  which  is  a  system  of  first-order  ordinary 
differential  equations.  One  approach  of  this  transformation 
which  is  mathematically  appealing  yet  easy  to  perform  is 
described  in  this  article. 
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Fig.  1.  Subdivision  of  the  ME  A  into  five  regions. 


2.  Dynamical  model  of  the  DMFC 

In  this  article  the  dynamical  behaviour  of  the  DMFC  is 
characterised  by  functions  of  place  z  and  time  t  for  methanol, 
oxygen,  carbon  dioxide  and  electrical  and  ionic  potential. 
Temperature  changes  and  water  flow  will  be  neglected,  since 
they  show  but  little  influence  on  the  dynamics  of  the  DMFC 
(see  [7]). 

The  model  of  the  dynamical  behaviour  of  the  DMFC 
is  restricted  to  its  main  part,  the  membrane  electrode 
assembly  (ME  A).  Generally,  the  ME  A  can  be  divided  into 
five  regions,  namely  two  catalytic  converter  layers,  where 
the  oxidisation  and  the  reduction  takes  place,  two  diffusion 
layers,  which  feed  the  catalytic  converter  layers  with  oxygen 
and  methanol,  and  a  membrane.  They  are  shown  in  Fig.  1. 
The  partition  into  five  layers  is  caused  by  several  physical 
parameters  due  to  different  materials.  The  modelling  of  a 
DMFC  is  restricted  to  one  spatial  dimension.  The  approach 
may  well  be  extended  to  a  two-dimensional  model,  yet  the 
error  of  a  one-dimensional  model  is  less  than  10%  (see 
[3]). 

Based  on  conservation  laws,  Rosenthal  [7]  obtained  a  set 
of  PDEs  for  the  quantities  concentration  of  methanol  CMe>  of 
oxygen  co2  and  of  carbon  dioxide  cqo2  as  well as  partial  pres¬ 
sure  of  oxygen  /?o2and  carbon  dioxide  pco2  •  Additionally, 


The  model  of  the  ME  A  uses  liquid  and  gas  phases,  both 
reacting  independently,  whereas  the  void  ratio  is  supposed 
to  be  uniform  at  the  whole  ME  A.  To  simplify  matters,  these 
equations  are  confined  to  reduction  of  oxygen  and  oxidation 
of  methanol  while  neglecting  intermediate  reactions.  Instead, 
the  global  kinetics  of  the  ideal  equation  of  reaction  is  used. 


2.7.  Differential  equations 

According  to  the  diffusion  law  by  Fick,  the  concentration 
of  methanol  CMeis  described  in  all  regions  by  diffusion  equa¬ 
tions: 


9cMe(F  z) 
dt 


d2CMe(*,  Z) 
dz2 


+  /Me  (6  Z) 


for  t  e  (0,  Tfland  T  e  R.  The  constants  hi  for  any  index  i  here 
and  in  all  following  equations  are  placeholders  for  physical 
parameters  of  the  DMFC.  They  may  vary  between  different 
layers.  A  complete  account  of  equations  and  constants  can 
be  found  in  [7]. 

Inevitably,  no  known  membrane  inhibits  methanol 
crossover  completely.  Hence,  methanol  passes  through  the 
membrane  by  diffusion,  convection  (the  whole  fall  of  pres¬ 
sure  difference  takes  place  in  the  membrane)  and  migration 
of  the  protons  (molecules  of  methanol  are  dragged  through 
the  membrane),  and  thus  may  interchange  between  anodic 
and  cathodic  layers.  These  additional  effects  are  reflected  by 
three  terms  in  /CMe  in  the  region  73 .  The  oxidation  of  methanol 
takes  place  in  the  catalytic  converter  layers  /2and  I4.  Hence, 
the  right-hand  side  of  Eq.  (1)  is  given  by 


(  0 


/Me  =  < 


k2  exp(k3(<pa  -<Pm~  h))CMe 

9  (j  dco2  dcco2  ,  7  d(Pm 

—  CMe  — - b  &20  — - b  fei^— 

OZ  \  oz  OZ  OZ 

k5  exp(k3(<pk  -<pm-  k4))cMe 

0 


z  e  I 2 

z^zeh 

z  e  I 4 
<=  z  e  1 5 . 


the  model  includes  boundary  value  problems  of  the  electrical 
potentials  <paand  cp^as  well  as  of  the  ionic  potential  cpm.  The 
spatial  domains  of  the  differential  equations  are  the  regions 
shown  in  Fig.  1 .  Their  physical  dimensions  for  the  MEA  be¬ 
longing  to  the  DMFC  described  in  [7]  are 

h  =  (zi,Z2)  =  (0,  520  p,  m), 
h  =  (Z2,  zs)  =  (520,  550  \x  m), 
h  =  (Z3,  za)  =  (550,  760  \x  m), 

14  =  (Z4,  Z5)  =  (760,  810  [x  m), 

15  =  (Z5,  Z6)  =  (810,  1330  \x  m). 


Inside  the  membrane  oxygen  and  carbon  dioxide  are  de¬ 
scribed  by  concentrations  co2and  cco2-  The  PDEs  for  z  £ 
/3 are  similar  to  Eq.  (1)  with  /o2  =  fco2  =  0-  Outside 
the  membrane  both  gases  are  described  by  partial  pres¬ 
sures  /?o2and  Pco2  •  Their  dynamics  is  represented  with 
t  £  (0,  Tfland  z  £  h  U  h  U  I 4  U  75 by  four  nonlinear  PDEs: 


9po2,co2(6  z) 
St 


9 

—  N 
dz 


(po2,  Pco2 


+  /o2,co2(6  z). 


The  nonlinear  transport  term  N(po2 ,  pco2 ,  dpo2/dz,  dpco2  / 
dz) in  Eq.  (3)  describes  the  transport  of  oxygen  and  carbon 
dioxide  outside  the  membrane  by  a  Binary-Friction-model. 
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The  nonlinear  terms  /o2and  /co?are  given  by 


fOo  =  < 


ro 

k22  exp(£8(<pa 

k22  exp(&8(</Jk 

0 


fm 


k9))Po7 

k9))PO 


and 


fco2  =  < 


ro 

k24  exp(&3(<pa 

k25  exp(£3(<pk 

0 


fm 


Z  €  I\ 
Z  €  I 2 

z  e  I 4 
ze  I 5 


h))CMe 

h))CMe 


z  e  h 
ze  I 2 
ze  I 4 
z  e  /5, 


(4) 


(5) 


which  represent  in  /2and  /4the  reduction  of  oxygen  and  the 
synthesis  of  carbon  dioxide  by  oxidation  of  methanol,  respec¬ 
tively. 

The  electrical  and  ionic  potentials  ^can  be  described  sim¬ 
ilarly  to  CMe-  Yet  the  dynamics  of  (pis  much  faster  than  the 
changes  of  CMe-  The  dynamics  of  the  potentials  may  therefore 
be  neglected,  so  that  their  behaviour  is  described  solely  by 
boundary  value  problems.  Thus,  the  potentials  (pare  described 
by  time-independent  PDEs: 

o  =  ki6-  J  ^  +  /ya(z),  z€/iU/2,  (6) 

d2cpm(z) 

0  =  *17  ^2  +  U mfe)’  z  €  h  u  h  U  h  and  (7) 

d2(Pk(z) 

0  =  *18^^  +  /**(*),  zel 4  U  /5.  (8) 

oz 

The  nonlinear  terms: 


are  necessary  for  every  PDE  depending  on  time  t  in  each  re¬ 
gion.  Additionally,  every  time-independent  PDE  needs  two 
spatial  boundary  conditions  for  each  region. 

The  initial  values  for  the  PDEs  in  Section  2.1  can  be  cho¬ 
sen  arbitrarily.  For  instance,  po2(0,  z) and  pco2( 0,  z)as  initial 
values  for  Eq.  (3)  reflect  merely  partial  pressures  at  t  =  0. 
They  are  determined  by  certain  environmental  conditions. 

The  boundary  conditions  may  be  separated  into  two 
groups.  The  boundary  values  of  the  MEA  at  z  =  z\  =0  p, 
mand  z  =  Z6  =  1330  p,  mare  determined  again  by  some  en¬ 
vironmental  conditions  and  may  be  chosen  arbitrarily.  In 
contrast,  the  inner  conditions  at  the  layer  contacts  at  z  = 
z 2,  . . . ,  zshave  to  reflect  the  physical  properties  of  the  MEA. 
In  general,  any  solution  of  the  PDEs  must  be  continuous  on 
the  entire  MEA.  Hence,  the  inner  boundary  conditions  are 
equalities  between  left-hand  and  right-hand  sides.  For  po2in 
Eq.  (3)  for  instance,  the  inner  boundary  conditions  read: 

P02,\(t,Zi)  =  P02,r(*>  Zi),  2=2,...,  5  (12) 

for  t  e  (0,  T].  Here  as  well  as  in  the  entire  sequel,  the  index 
lmarks  the  left-hand  value  at  zzand  raccordingliy  the  right- 
hand  value  at  Zi : 

Po2,i(*>  zd  =  lim  po2(t,  Zi  -  s )  and 

£'—>() 

Po2Ak  Zi)  =  lim  po2(t,  zi  +  e),  /  =  1 . 5. 

£ — >  0 

Additionally,  the  right-hand  sides  of  the  PDEs  of  Section 
2.1  describing  the  material  flow  must  be  at  least  continuous. 
Therefore,  additional  inner  boundary  conditions  are  given  by, 
here  exemplarily  in  case  of  Eq.  (3): 


f(Pa 


f fim 


0  Z  ell 

k6  exp(k3(^a  -<Pm-  &4))cMe  +  k7  (<P&  ~  <Pm  ~  k9))P02  Z  e  h, 

'  ho  exp(k3(^a  -  (Pm  ~  &4))cMe  +  *11  exp(kg(^a  ~  <Pm  ~  k9))P02  quadz  €  h 
0  quadz  e  /3 

I  k\2  exp (k3(cpk  -<pm-  k4))cMe  +  h3  exp(k8(^k  -<Pm-  k9))po2quadz  e  /4, 


(9) 

(10) 


k  14  exp(k3(<pk  -  cpm  -  k4))cMe  +  k\5  Qxp(h(cpk  -  (pm  -  k9))po2  ze  u 
0  z  e  Is 


(11) 


reflect  the  potentials  which  depend  on  CMeand  po2. 

All  in  all,  the  dynamic  model  of  the  DMFC  is  given  by 
three  time-dependent  PDEs  for  each  region  I\ ,  . . . ,  Is,  a  sin¬ 
gle  time-independent  PDE  for  each  of  the  regions  I\,  /3and 
^and  two  time-independent  PDEs  for  each  of  the  regions 
/2and  /4.  Ten  of  the  time-dependent  PDEs  are  nonlinear.  Un¬ 
fortunately,  each  quantity  has  a  nonlinear  right-hand  side/in 
at  least  one  region. 

2.2.  Boundary  and  initial  conditions 

To  characterise  the  dynamic  behaviour  of  the  DMFC  com¬ 
pletely,  two  initial  conditions  and  two  boundary  conditions 


N\  —  Nr  at  z  =  zi,  •  •  • ,  £5- 

Generally,  with  y  being  a  placeholder  for  concentrations  and 
partial  pressures,  the  set  of  conditions  can  be  summarised  as 

Bv y(t,  Zi)  =  yT(t)  and  B\y(t,  Zi+i)  =  y\(t)  (13) 

with  operators  £iand  Br for  each  differential  equation  and 
/  =  1,  . . . ,  5.  The  time-independent  PDEs  describing  the 
electrical  and  ionic  potentials  have  inner  boundary  conditions 
similar  to  Eq.  (13),  except  that  they  are  time-independent, 
too. 
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3.  State  space  representation  of  the  DMFC 


The  dynamics  of  the  DMFC  is  represented  by  15  nonlinear 
and  7  linear  partial  differential  equations  of  degree  two.  This 
model  of  the  DMFC  is  not  well  suited  in  a  control-theoretical 
sense,  i.e.  to  analyse  the  cells  dynamic  behaviour  and  to  syn¬ 
thesise  a  proper  controller.  Therefore,  these  differential  equa¬ 
tions  will  be  transformed  into  an  equivalent  mathematical  for¬ 
mulation,  the  state  space  representation.  A  state  space  model 
allows  detailed  analysis  of  the  cells  temporal  and  spatial  be¬ 
haviour,  optimisation  of  physical  parameters  of  the  cell  and 
is  well  suited  for  synthesis  of  a  controller  to,  e.g.  stabilise  the 
cells  voltage  or  reject  environmental  disturbances.  Besides, 
there  are  powerful  tools  for  simulation,  analysis,  synthesis 
and  optimisation  of  the  dynamic  behaviour  of  a  state  space 
model. 

The  transformation  of  the  differential  equations  devel¬ 
oped  in  Sections  2.1  and  2.2  into  a  state  space  model,  i.e.  a 
set  of  linear  ordinary  first-order  differential  equations  will  be 
carried  out  in  several  steps.  Laplace  transforms  are  applied 
to  eliminate  partial  time  derivatives.  The  result  is  a  set  of 
ordinary  differential  equations  either  of  place  z  or  of  both 
frequency  s  and  place  z.  Afterwards,  the  dependency  upon 
place  z  is  eliminated  by  the  method  of  weighted  residuals, 
that  is  by  approximating  the  spatial  dependency  as  accurately 
as  required.  Finally,  the  nonlinearities  of  the  differential 
equations  are  treated  to  yield  a  numerically  convincing  yet 
linear  approximation  of  the  dynamics  of  the  DMFC,  which 
is  the  state  space  model,  i.e.  a  system  of  ordinary  linear 
second-order  initial  value  problems  [5]. 

3.1.  Laplace  Transform 


With  the  Laplace  transform: 


*oo 


Y(s,z)  =  Cy(t,  z)  =  I  y(t,z)c  51  dt, 


(14) 


Jt= o 

linear  partial  differential  equations  may  be  transformed  into 
a  boundary  value  problem: 

DY(s,  z )  =  Upis,  z),  z  €  f  =  (zi,  zi+i)  C  M, 


(15) 


with  D  being  a  linear  operator,  which  depends  on  the  complex 
frequency  parameter  s  (see  [4]).  Capital  letters  denote  the  cor¬ 
responding  time-dependent  lowercase  letters  in  the  frequency 
domain,  a  notation  which  is  adopted  in  the  entire  sequel.  Ap¬ 
plying  the  Laplace  transform  (14)  to,  for  instance,  Eq.  (1), 
the  partial  differential  equation  in  CMe(6  z)c an  be  written  as 
an  ordinary  boundary  value  problem  in  Cm&(s,  z) and  reads: 


sCMq(s,  z) 


d2CMe(s,  Z) 

3  z2 


£/Me(*,  Z)  +  CMe(0,  z). 


(16) 


The  right-hand  side  of  Eq.  (16)  is  the  equivalent  of  U /.in  Eq. 
(15).  U contains  all  possibilities  to  exert  influence  within  the 


region  i  =  1, . . . ,  5.  A  description  similar  to  (15)  can  be 
stated  for  all  time-dependent  PDEs  introduced  in  Section  2.1. 

A  precondition  for  applying  the  Laplace  transform  to  a 
PDE  is,  that  neither  mixed  partial  derivative  occurs.  For  most 
physical  systems  this  precondition  is  satisfied.  Here  too,  all 
PDEs  introduced  in  Section  2.1  include  either  spatial  or  tem¬ 
poral  partial  derivatives,  but  not  mixed  ones. 


3.2.  Method  of  weighted  residuals  and  collocation 

One  of  several  possibilities  to  transform  partial  differen¬ 
tial  equations  or  boundary  value  problems  is  the  method  of 
weighted  residuals  [4] .  The  fundamental  idea  of  this  method 
is  the  approximation  of  the  exact  solution  y  by  a  linear 
combination: 

N  N 

y(t,  z)  =  T,yj(t)<Pj(z)  or  F(s,  z)  =  Yj(s)<pj(z ) 

7=1  7=1 

(17) 


of  N  e  Nbasis  functions  <pj(z).  Thus,  solving  PDEs  or 
boundary  value  problems  in  y(t ,  z)and  Y(s ,  z)is  replaced  by 
calculation  of  the  coefficients  yj(t) and  Yj(s),  respectively. 
The  basis  functions  <pj(z) in  Eq.  (17)  can  be  arbitrarily 
chosen,  they  merely  need  to  be  linearly  independent. 

For  modelling  the  DMFC,  Chebyshev  polynomials  have 
been  used  as  basis  functions  0;-(z).  They  are  an  orthogonal  set 
of  functions,  which  will  play  a  major  role  later  (see  Section 
3.3).  Furthermore,  approximations  based  on  Chebyshev  poly¬ 
nomials  converge  rather  fast,  hence  only  few  basis  functions 
are  necessary  and  N  remains  small.  Not  at  least,  calculations 
based  on  simple  polynomials  are  computationally  cheap  [6]. 

The  main  advantage  of  the  method  of  weighted  residuals  is 
a  small  number  of  basis  functions  A  compared  to  several  other 
approaches.  This  leads  to  a  comparatively  small  dimension  of 
the  resulting  system  of  linear  equations.  Moreover,  the  solu¬ 
tion  is  approximated  not  only  at  certain  discretisation  points. 
The  coefficients  }y(0represent  y(t,  z) at  every  point  in  the  re¬ 
gion  likewise  for  Laplace  transformed  quantities.  Thus, 
interpolation  is  avoided  when  evaluating  y(t,  z)or  Y(s ,  z)at 
some  z  e  /* ,  i  =  1, . . . ,  5,  which  is  particularly  important  for 
system  analysis. 

The  coefficients  yj{t) and  F/(^)have  to  be  calculated  re¬ 
specting  two  constraints.  Namely,  the  approximating  func- 
tion  yor  yresp.  has  to  satisfy  the  boundary  value  problem 
(15)  as  well  as  the  corresponding  boundary  conditions  (13). 

The  coefficients  yj(t) of  Eq.  (15)  are  calculated 
by  minimisation  of  the  approximation  errors  Rn  = 
Dy(t,  z)  —  Dy(t,  z) and  R n  =  DY(s,  z)  —  DY(s,  z)resp.,  i.e. 
the  weighted  residuals: 

r?d+\ 

rjc=  RNWk(z)dz,  k  =  1, . . . ,  N  -  2,  i  =  1, ...  5, 

Jzi 


(18) 
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should  vanish.  With  a  set  of  N  —  2weighting  functions 
Wk(z ),  (18)  is  a  system  of  N  —  2equations  to  determine  the 
N  coefficients  yj(t) or  Yj(s).  Two  equations  are  yet  required 
to  uniquely  solve  (18).  They  are  given  by  the  boundary 
conditions  (see  Eqs.  (21)  and  (22)  below). 

The  weighting  functions  Wk(z) can  be  arbitrarily  cho¬ 
sen,  they  only  need  to  be  linearly  independent.  The  sim¬ 
plest  choice  is  the  Dirac  impulse  at  some  point  zCp,k,  which 
is  called  collocation,  pseudospectral  or  method  of  selected 
points,  Wk(z)  =  8(z,  zCp,k)fov  k  =  1, . . . ,  N  —  2.  Now  the 
weighted  residuals  (18)  reduce  to  r f  =  Rn(s,  zCp,k)ov  r f  = 
Rn(U  Zcp,yt)for  k  =  1,  . . . ,  N  -  2. 

For  boundary  value  problems  without  time-dependency 
follows: 

N 

D</>j(zcp,k)  =  uiiizcp'k),  k  =  1, . . . ,  N  —  2.  (19) 

^  =vkj  =Uk 

For  ^-depending  boundary  value  problems  the  method  of 
weighted  residuals  for  k  =  1 ,  . . . ,  N  —  2results  in  linear 
equation  systems  depending  on  s: 


Eq.  (19)  resp.  (20)  together  with  Eq.  (21)  resp. 
(22)  completely  characterise  the  approximation  y(z)resp. 
Y(s,  z)according  to  Eq.  (17)  of  the  solution  of  the  DMFCs 
model  introduced  in  Section  2.1. 

3.3.  Linear  state  space  model  of  the  DMFC 

For  all  seven  time-independent  boundary  value  problems 
we  recall  Eq.  (19)  and  define  the  quantities: 

y  ~  ( yj )  G  R1N ,  u  =  (ui)  G  and 

V  =  (vkj)  e  R7NxlN. 

On  the  other  hand,  for  each  of  the  three  frequency-depending 
boundary  value  problems,  we  define  according  to  Eq.  (20): 

Y(s)  =  ( Yj(s ))  e  R5N,  U(s )  =  (Ui(s))  e  R5N  and 

V  =  (vkj)  e  R5Nx5N,  W  =  (■ wtj )  €  R5Nx5N. 

Now,  time-independent  boundary  value  problems  (19)  can 
simply  be  written  as  a  linear  matrix  equation: 


N 


E  w 

j= i 


\ 


ks  f  j(zCp,k)  T  DZ<P j(Zcp,k) 


—  U ifs,  Zcp,k )  • 

s - v - / 

=Uk(s ) 


Vy  =  u.  (23) 

Similarly,  frequency-depending  boundary  value  problems 
(20)  now  simply  read: 


(2°)  (ksW  +  V)f  (s)  =  U(s). 


(24) 


Here,  the  operator  D  is  replaced  by  a  place-depending,  time- 
invariant  operator  Dz  and  a  time- variant  operator.  Since  all 
PDEs  of  the  DMFC  are  of  degree  one  for  the  derivation  with 
respect  to  time,  the  time- variant  operator  can  be  written  as 
Csand  k(d/dt ),  respectively. 

/V 

With  yor  Tresp.  satisfying  the  boundary  value  problem 
(15),  they  still  have  to  fulfill  the  corresponding  boundary 
conditions  (13).  According  to  Section  2.2,  all  inner  and  outer 
boundary  conditions  at  z  =  z\,  . . . ,  Z6can  be  written  in  a  stan¬ 
dardised  form  with  operators.  They  shall  be  transformed  sim¬ 
ilarly  to  Eq.  (19)  in  case  of  a  boundary  value  problem  and 
Eq.  (20)  for  a  partial  differential  equation. 

Applying  Eq.  (17)  to  Eq.  (13)  yields: 

N  N 

B r  ^2  yj<t>j(zi)  =  y r  and  B\  ^  yjfj(zi+\)  =  y\, 

7=1  7=1 

i=  1,...,5  (21) 

for  time-independent  boundary  value  problems  and 

N 

Br  Y,  Yj(s)<pj(z.i  )  =  Yr(s)  and 
7=1 

N 

B\  Y j{s)(j> j{zi+\)  =  i)(s),  i  =  1 . 5,  (22) 

7=1 


Eqs.  (23)  and  (24)  are  a  linear  state  space  model  of  the  DMFC. 
Yet,  for  efficient  analysis,  optimisation  and  synthesis  of  the 
dynamical  system  further  transformations  are  necessary  to 
yield  a  simple  input/output  relation. 

The  rows  of  V,  Eq.  (23),  are  created  either  by  the  method 
of  weighted  residuals  (in  total  7 (N  —  2)rows)  or  by  boundary 
conditions  (7  x  2rows).  The  latter  ones  are  pairwise  linearly 
independent,  if  only  the  boundary  conditions  at  ziand  Z6&re 
of  Dirichlet  type.  This  applies  to  the  DMFC  as  can  be  seen 
in  Eq.  (21).  The  1(N  —  2)rows  of  V created  by  the  method  of 
weighted  residuals  are  pairwise  linearly  independent,  if  the 
set  of  basis  functions  <pj(z ),  j  =  1,  . . . ,  A^(see  Eq.  (17)),  is 
orthogonal.  Since  both  sets  of  rows  are  necessarily  linearly 
independent,  the  matrix  V is  regular  for  orthogonal  basis  func¬ 
tions. 

With  a  regular  matrix  V,  Eq.  (23)  reads: 
y=V~1u 

with  ybeing  the  unique  solution  for  the  corresponding  time- 
independent  boundary  value  problems. 

Frequency-depending  boundary  value  problems  cannot  be 
transformed  similarly,  since  the  left-hand  side  of  Eq.  (24) 
depends  on  the  frequency  parameter  s.  Though,  Eq.  (24)  may 
be  further  transformed  to 

1  , 

sY(s)  =  -  W~\-VY(s )  +  U(s)) 
k 


for  frequency-depending  boundary  value  problems. 


(25) 
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Fig.  2.  Block  diagram  of  the  state  space  model  of  linear  partial  differential 
equations. 

if  only  W  is  regular.  For  general  boundary  conditions,  and 
boundary  conditions  of  the  DMFCs  model,  too,  the  regular¬ 
ity  of  W  cannot  be  guaranteed.  Nevertheless,  the  boundary 
conditions  of  the  DMFC  may  be  transformed  to  a  differential 
equation,  if  they  are  differentiated  in  time,  and  thus  create  a 
regular  matrix  W. 

Reverse  Laplace  transform  of  Eq.  (25)  yields  a  linear  state 
space  model  of  time-dependent  boundary  value  problems: 

}’(t)  =  \w~'(-Vy(t)  +  fi(0),  t  e  (0,  T],  (26) 

k 

Altogether,  the  dynamical  behaviour  of  the  DMFC  is  captured 
by  a  state  space  model  (26)  with 

1 

system  matrix  :  A  =  —  W  V  and 

k 

1  , 

input  matrix  :  B  =  -W  1 .  (27) 

k 

Fig.  2  shows  the  corresponding  block  diagram  for  Eq.  (26). 
The  input  vector  uis  given  by  5 (N  —  2) values  uft) included  in 
U(s,  z)at  collocation  points  zcp,&(see  Eq.  (15)).  The  remain¬ 
ing  5  x  2 values  are  the  boundary  conditions,  i.e.  the  deriva¬ 
tive  of  the  right-hand  sides  of  Eq.  (13).  The  output  vector 
y(f)consists  of  coefficients  yj(t),  which  in  fact  represent  the 
solution  of  the  partial  differential  equation. 

3.4.  Linear  time-variant  state  space  model  of  the  DMFC 

The  results  of  Section  3.3  do  not  yet  incorporate  any  non¬ 
linearity  of  the  partial  differential  equations  introduced  in 
Section  2.1.  To  yield  a  linear  state  space  model  for  analysis, 
optimisation  and  control  of  the  DMFC,  all  nonlinear  terms 
must  be  eliminated.  Of  course,  any  kind  of  linearisation  of 
the  differential  equations  is  straightforward.  But  a  priori  no 
guaranteed  error  bounds  can  be  stated  explicitely.  To  reduce 
the  linearisation  errors  and  particularly  keep  the  solution  in 
tight  error  bounds,  linear  approximations  of  the  nonlinear 
terms  are  calculated  iteratively. 

Here  again,  the  method  of  weighted  residuals  is  proposed. 
Two  possibilities  are  at  hand.  The  first  one  is  a  classical  New¬ 
ton  iteration,  in  which  the  nonlinear  problem  is  discretised 
with  the  method  of  weighted  residuals.  The  resulting  nonlin¬ 
ear  matrix  equation  can  be  solved  iteratively  by  linearisation 
(see  e.g.  [9]).  The  second  possibility  is  to  linearise  about  a 
function  and  discretise  the  linear  problem  afterwards.  The 
second  variant  is  the  Newton-Kantorovich  method.  It  is  pre¬ 


ferred  in  this  context  [2],  since  it  is  computationally  much 
cheaper  and  numerically  more  stable. 

Yet  there  is  another  possibility  to  yield  a  state  space  model 
in  a  control-theoretical  framework.  Accepting  a  slight  change 
in  the  notion  of  states,  the  nonlinearities  in  partial  differen¬ 
tial  equations  similar  to  Eq.  (1)  or  (3)  may  be  considered  as 
time-depending  variation  of  the  states,  i.e.  as  auxiliary  inputs 
ua(t,  z).  Applying  Laplace  transform  and  collocation,  they 
shall  read  Ua(s)  =  Cu^t,  zCp,k)at  collocation  points  zCp,k, 
k  =  1, . . . ,  N  —  2.  Of  course,  no  external  influence  on  the 
DMFC  can  be  exerted  with  these  auxiliary  inputs.  Instead, 
all  additional  inputs  are  already  determined  by  the  partial 
differential  equations. 

With  these  auxiliary  inputs,  Eq.  (25)  will  be  extended  to 

sY(s)  =  lw-\-VY(s )  +  U(s)  +  U.JS)).  (28) 

k 

The  matrices  V,  W and  the  input  term  fiaare  the  same  as  in  (25). 
The  inputs  Ua(s) are  completely  determined  by  the  PDEs. 
Consequently,  their  dependency  upon  time  or  frequency  resp. 
can  be  rewritten  as  time- variant  dependency  upon  the  states 
of  the  state  space  model  of  the  DMFC.  As  an  example,  the 
concentration  of  methanol  CMe(see  Eq.  (1)),  shall  be  consid¬ 
ered  in  the  anodic  catalytic  converter  layer  /2(see  Fig.  1  in 
Section  2).  There,  fia(0reduces  to  a  term  14(0  ^depending  on 
time  t.  The  entries  of  Va(t)a re  given  by 

( yf)k,j  —  kl  exp(k3(^afecp,k)  —  ^m(^cp,k)  —  ^4))0jfecp,k) 

for  k  =  1 ,  . . . ,  N  —  2and  j  =  1 ,  . . . ,  N.  The  entries  still 
missing,  i.  e.  the  remaining  two  rows  of  14(0,  are  determined 
by  the  corresponding  boundary  conditions  (see  Section  3.3). 

Similarly  to  Eq.  (25),  the  reverse  Laplace  transform  of 
Eq.  (28)  yields  a  linear  state  space  model  of  the  nonlinear 
boundary  value  problems: 

y(t)  =  \w~\-V%t)  +  Va(t)y(t)  +  fl(0),  t  e  (0,  T], 
k 

(29) 

The  dynamical  behaviour  of  the  DMFC  represented  by  a  state 
space  model  (29)  has  an  input  matrix  similar  to  Eq.  (27). 
The  system  matrix  is  divided  into  a  time-independent  part 
W~l  Land  a  time-depending  part  W~l  Va(0,  i.e. 

1  -i  -i 

system  matrix  A  =  -(— W  V  -\-W  va(0)  and 

,  k  (30) 

input  matrix  B  =  -W~l . 

k 

The  corresponding  block  diagram  for  Eq.  (27)  in  conjunc¬ 
tion  with  Eq.  (28)  can  be  seen  in  Fig.  3. 

4.  Results  and  discussion 

The  state  space  model  of  the  entire  DMFC  is  a  compo¬ 
sition  of  system  and  input  matrices.  Each  partial  differential 
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Fig.  3.  Block  diagram  of  the  state  space  model  of  nonlinear  partial  differ¬ 
ential  equations. 

equation  provides  a  state  space  model  according  to  Eq.  (26)  or 
(29)  resp.  Consequently,  the  block  diagram  of  the  state  space 
model  of  the  DMFC  is  composed  of  several  block  diagrams 
shown  in  Figs.  2  and  3. 

To  validate  the  approach  of  modelling  the  DMFC  in  a 
control-theoretical  framework,  the  state  space  model  has  been 
simulated  for  several  environmental  conditions,  i.e.  bound¬ 
ary  and  initial  values.  The  results  have  been  compared  with 
simulation  results  of  the  FEM  model  of  the  PDEs  developed 
by  Rosenthal  [7].  Additionally,  measuring  data  of  the  DMFC 
described  by  Rosenthal  have  been  taken  into  account  if  avail¬ 
able. 

An  important  and  widely  used  quality  of  a  fuel  cell  is  its 
current  density-voltage  characteristic.  It  has  been  simulated 
and  compared  for  various  environmental  conditions.  Exam- 
plarily,  a  state  space  simulation,  a  FEM  simulation  and  mea¬ 
suring  data  for  a  single  set  of  boundary  conditions,  which 
are  specified  in  the  caption,  are  presented  in  Fig.  4.  As  can 
be  seen,  the  state  space  model  closely  resembles  the  FEM 
model.  Even  though  in  the  state  space  model  the  temperature 
is  fixed,  the  differences  between  the  FEM  model  and  the  state 
space  model  are  negligible.  The  errors  between  the  models 
and  the  original  data  are  caused  be  the  simplifications  listed 
in  Section  2. 

The  current  density-voltage-characteristic  plots  the  cells 
voltage  against  the  steady  state  working  point  of  the  current 
density  for  constant  boundary  values  and  physical  param¬ 
eters.  It  allows  the  analysis  of  effects  certain  parameters 
show  upon  the  performance  of  the  DMFC  when  operated 
stationarily.  Subsequently,  environmental  conditions  as  well 
as  physical  parameters  can  be  altered  to  optimise  particular 
properties  of  the  cell. 

Yet  the  state  space  model  is  particularly  well  suited  for 
dynamical  simulations.  The  step  response  is  a  standard 
method  to  analyse  the  systems  dynamical  behaviour.  An 
example  is  shown  in  Fig.  5.  The  methanol  concentration 
Qvie(£i)at  the  anodic  outer  boundary  (solid  step  function)  is 
stepwise  altered  at  time  t  =  Oand  60  s.  The  plot  shows  the 
reaction  of  the  state  space  model  (solid  line)  and  the  reaction 
of  the  FEM  model  (dashed  line)  on  this  varying  environ¬ 
mental  condition.  Here  again,  the  state  space  model  closely 


Fig.  4.  current  density-voltage  characteristic;  applied  boundary  conditions 
are  pco2(zi)  =  1  bar,  pco2{ze)  =  Obar,  po2(z\)  =  Obar,  po2(ze)  =  Ibar 
and  CMeCi)  =  750 mol m~3,  CMe(^6)  =  0molm~3. 

resembles  the  FEM  model.  Which  of  both  models  reflects 
the  cells  dynamics  more  precisely  cannot  be  evaluated,  since 
no  measuring  data  are  available  for  pure  steps  of  methanol 
concentration. 

Yet  the  main  advantage  of  the  state  space  model  cannot 
be  shown  in  figures.  Primarily,  the  state  space  model  as  a 
standard  approach  of  control  theory  allows  the  analysis  and 
optimisation  of  the  cells  dynamics.  For  instance,  improving 
the  behaviour  of  the  cell  in  case  of  rapidly  varying  electri¬ 
cal  load  will  be  part  of  future  work.  The  behaviour  of  the 
cell  may  also  be  improved  by  a  closed  loop  controller,  e.g.  to 
stabilise  the  cells  voltage  by  adaptation  of  methanol  concen¬ 
tration.  Again,  a  state  space  model  is  the  preferred  description 
a  controller  synthesis  is  based  upon. 

Another  advantage  of  the  state  space  model  is  the 
possibility  of  choosing  the  boundary  values  freely.  The  FEM 
model  proposed  by  Rosenthal  [7]  is  restricted  to  certain  com¬ 
binations  of  boundary  values  due  to  its  inner  structure.  The 


Fig.  5.  Reaction  of  the  cell  to  an  anodic  step  of  methanol  concentration  with 
boundary  conditions  pco2(t,  zi)  =  1  bar,  pco2(t,  zf)  =  Obar,  po2(t,  zi)  = 
Obar,  po2(t,  Z6 )  =  1  bar,  CMefi,  £6)  =  0molm~3  and  cp^t,  zf)  =  0  V. 
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inverse  current  density-voltage  characteristic,  for  instance, 
can  be  calculated  directly  only  by  the  state  space  model. 

5.  Conclusions 

In  this  article,  an  approach  is  shown,  which  transforms 
the  dynamical  model  of  a  direct  methanol  fuel  cell  into  a 
state  space  representation.  This  representation  is  particularly 
well  suited  to  analyse  the  systems  dynamical  behaviour,  to 
optimise  certain  of  its  physical  properties  and  to  synthesise  a 
controller  when  necessary.  Therefore,  the  state  space  model  is 
an  important  tool  to  develop  robust  and  efficient  applications 
of  fuel  cells.  For  instance,  synthesis  of  a  controller  to  ensure 
stable,  robust  and  efficient  operation  of  a  direct  methanol  fuel 
cell  will  be  part  of  future  work. 
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